echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 4.1 (Evaluation of a Taylor Series).
% Section	4.1,	Taylor Series and Calculation of Functions, Page 203
echo on; clc; format long; hold off; clear

% This program investigatges Taylor approximations.

%      Pn(x) = c(1) + c(2)x + c(2)x^2 + ... + c(n+1)x^n

% where the degree n of approximation is large  (n ~ 25).

% Coefficient lists for several functions have been

% stored in M-files named;   zcos  zsin  ztan  zexp

% zacos  zasin  zatan  zcosh   zsinh   zsqrt   zlog
 
% zsqrt4   zinv   zemx2d2   zcosde   zsinde   zlogq

pause	% Press any key continue.

clc;
%       Approximations for  tan(x)

% Issue the command  ztan  to load the coefficients

% into the array  C.  The function name is loaded

% into the variable fun, the degree is loaded into  m.

% The endpoints of [a,b] are loaded into  a  and b.

% Load the Taylor coefficients.

[fun,dfun,ifun,x0,m,C,Ax] = ztan;

pause 	% Press any key continue.
clc;
a = Ax(1,1);    % You can change the left  endpoint a.
b = Ax(1,2);    % You can change the right endpoint b.
c = Ax(1,3);
d = Ax(1,4);
n = m;          % You can change the degree n.
C = flipud(C);
C = C(1:n+1);
C = flipud(C);

pause	% Press any key continue.

clc; clg;
h = (b-a)/200;
X = a:h:b;
x = X; 
Y = eval(fun);
axis([a b c d]);...
P = polyval(C,X);...
plot(X,Y,'-g',X,P,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = ['Comparison of ',fun,' and P'];...
Mx2 = [Mx1,num2str(n),'(x)'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause	% Press any key to continue.

Mx1 = 'The function is  f(x) = ';
Mx2 = 'The interval is  ';
Mx3 = 'Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
Mx4 = 'The degree is  n = ';
Mx5 = ', and the coefficients list  C  is:';
clc,format short e,echo off,diary output,...
disp(''),disp([Mx1,fun]),...
disp([Mx2,'[',num2str(a),' , ',num2str(b),']']),...
disp(Mx3),disp([Mx4,num2str(n),Mx5]),...
for i=1:5:n+1, disp(C([i:min(i+4,n+1)])'); end,...
diary off, echo on

pause	% Press any key to graph f(x) - Pn(x).

clc; clg;
a = -1;    % You can change the left  endpoint a.
b =  1;    % You can change the right endpoint b.
h = (b-a)/200;
X = a:h:b;
x = X; 
Y = eval(fun);
P = polyval(C,X);
c = min(Y-P);
d = max(Y-P);
axis([a b c d]);...
plot(X,Y-P,'-r');...
hold on;...
plot([a b],[0 0],'b',[0 0],[c d],'b');...
xlabel('x');...
ylabel('y');...
Mx1 = ['The error;  ',fun,' - P'];...
Mx2 = [Mx1,num2str(n),'(x)'];...
title(Mx2);...
grid;...
axis;...
hold off;...
shg; pause	% Press any key for a list of numerical computations.

clc; format long;
a = -1.2;
b =  1.2;
X = a:0.1:b;
x = X;
Y = eval(fun);
P = polyval(C,X);
points = [X;Y;P;Y-P];
Mx1=['Taylor polynomial approximation of f(x) = ',fun];
Mx2='     x(k)               f(x(k))            Pn(x(k))           error';
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
diary off,echo on
